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Abstract - We study the ground-state (gs) phase diagram of the frustrated spin-1 J\—Ji antifer- 
romagnet with J2 = kJi > ( Ji > 0) on the honeycomb lattice, using the coupled-cluster method. 
We present results for the ground-state energy, magnetic order parameter and plaquette valence- 
bond crystal (PVBC) susceptibility. We find a paramagnetic PVBC phase for k c1 < k < k C2 , 
where « C1 « 0.207 ±0.003 and « C2 ~ 0.385 ±0.010. The transition at n ci to the Neel phase seems 
to be a continuous deconfined transition (although we cannot exclude a very narrow intermediate 
phase in the range 0.21 < k < 0.24), while that at k C2 is of first-order type to another quasiclassi- 
cal antiferromagnetic phase that occurs in the classical version of the model only at the isolated 
and highly degenerate critical point n = j. The spiral phases that are present classically for all 
values k > h are absent for all k < 1. 
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Two-dimensional (2D) frustrated quantum spin-lattice 
systems have become of huge interest both theoretically 
and experimentally [TH2]. Attention has particularly fo- 
cussed on the rich panoply of (zero-temperature, T = 0) 
quantum phase transitions that they exhibit [3j|4]. With- 
out thermal fluctuations the transitions are driven solely 
by the interplay of quantum fluctuations and any frustra- 
tion due to inherent competition between the interactions. 
Such frustration can arise either dynamically or geomet- 
rically. A prototypical example of the former is the well- 
studied J1-J2 Heisenberg antiferromagnet (HAFM) on the 
bipartite square lattice (see, e.g., Refs. [SHI] and references 
cited therein), where nearest-neighbour (NN) spins inter- 
act via a Heisenberg interaction with strength parameter 
J 1 > 0, which competes with a Heisenberg interaction with 
strength parameter J2 > between next-nearest neigh- 
bour (NNN) pairs. Similar prototypical models exhibiting 
geometrical frustration are the pure NN HAFMs on the 
triangular [5] and kagome lattices [3]. For cither form 
of frustration special interest then centres on the possi- 
ble appearance of novel quantum ground-state (gs) phases 
without the long-range order (LRO) that typifies the clas- 
sical gs phases of the corresponding models taken in the 
limit s — V co of the spin quantum number s of the lattice 
spins. Examples include various valence-bond crystalline 



solid phases and spin-liquid phases. 

Quantum fluctuations tend to be largest for the small- 
est values of s, for lower dimensionality D of the lattice, 
and for the smallest coordination of the lattice. Thus, for 
spin-i models in D = 2, the honeycomb lattice plays a 
special role. Frustration is easily incorporated via com- 
peting NNN and maybe also next-next-nearest-neighbour 
(NNNN) bonds. Such models and their experimental re- 
alisations have been much studied in recent years |10H16j . 
Additional interest has also sprung from the recent syn- 
thesis of graphene monolayers and other magnetic materi- 
als with a honeycomb structure. Theoretical interest was 
spurred by the discovery of a spin-liquid phase in the ex- 
actly solvable Kitaev model [T7|, in which spin-i particles 
reside on a honeycomb lattice. Hubbard models on the 
honeycomb lattice may also describe many of the relevant 
physical properties of graphene. For example, evidence 
has been found [TH] that quantum fluctuations are suffi- 
ciently strong to establish an insulating spin-liquid phase 
between the nonmagnetic metallic phase and the antiferro- 
magnetic (AFM) Mott insulator phase, when the Coulomb 
repulsion parameter U becomes moderately strong. For 
large values of U the latter phase corresponds to the pure 
HAFM on the bipartite honeycomb lattice, whose gs phase 
exhibits Neel LRO. However, higher-order terms in the 
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(a) Neel 



(b) spiral 



(c) anti-Neel 



Fig. 1: (Color online) The J\~Jn model on the honeycomb 
lattice (with J\ = 1), showing (a) the Neel, (b) spiral, and 
(c) anti-Neel states. The arrows represent spins located 
on lattice sites •. 



t/U expansion of the Hubbard model may lead to frustrat- 
ing exchange couplings in the corresponding spin-lattice 
limiting model, in which the HAFM with NN exchange 
couplings is the leading term in the large-C/ expansion. 
There is a growing consensus [TUl[TTl[r5rHB] that the frus- 
trated spin-i HAFM on the honeycomb lattice undergoes 
a frustration-induced quantum phase transition to a para- 
magnetic phase showing no magnetic LRO. 

Indirect experimental backup for such theoretical find- 
ings comes from recent observations of spin-liquid-likc 
behaviour in the layered compound BiaMn-tO^NOs) 
(BMNO) at temperatures below its Curie- Weiss tempera- 
ture [12]. In BMNO the Mn 4+ ions reside on the sites of 
(weakly-coupled) honeycomb lattices, but they have spin 
quantum number s = | . The successful replacement of the 
Mn 4+ ions in BMNO by V 4+ ions could lead to the real- 
ization of a corresponding s = h model on the honeycomb 
lattice. Other recent realizations of HAFMs on a hon- 
eycomb lattice include compounds such as NaaC^SbOg 
Ql] and InCu 2/3 V 1/3 3 [IS], in both of which the Cu 2 + 
ions in the copper oxide layers form a spin- ^ HAFM on a 
(distorted) honeycomb lattice. Others include the family 
of compounds BaM 2 (X0 4 )2 (M=Co, Ni; X=P, As) [21], 
in which the magnetic ions M are disposed in weakly- 
coupled layers where they reside on the sites of a hon- 
eycomb lattice. The Co ions have spins s = \ and the Ni 
ions have spins s = 1. Finally, recent calculations of the 
low-dimensional material /3-CU2V2O7 [22] show that its 
magnetic properties can be described in terms of a spin- \ 
model on a distorted honeycomb lattice. 

In recent papers [221121] we have studied the frustrated 
spin-i J1-J2-J3 model on the honeycomb lattice [TUlfTTl 
rT3"HTS] . Its Hamiltonian is given by 
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where index i runs over all honeycomb lattice sites, and 
indices j, k, and I run over all NN, NNN and NNNN sites 
to i, respectively, counting each bond once and once only. 
Each lattice site i carries a particle with spin s = | and a 
spin operator Sj = (sf,s^,sf). The lattice and exchange 
bonds are illustrated in Fig. [T] In earlier work we re- 
stricted ourselves to the case where J3 = J2 > 0, but we 



dealt both with the AFM case with Ji > [23] and the 
ferromagnetic (FM) case with J\ < [21] • Here we put 
J3 = 0, and hence consider the J1-J2 model. We restrict 
ourselves to the case where both bonds are of AFM type, 
J 1 > and J2 = kJi > 0. We henceforth set J\ = 1. 

The classical (s — > 00) gs phase diagram of the J\—J%- 
J3 model on the honeycomb lattice [TT][25] comprises six 
different phases when J\ > and the other two bonds, 
J2 and J3, can take either sign. Three are collinear AFM 
phases, one is the FM phase, and the other two are differ- 
ent helical phases (and see, e.g., Fig. 2 of Ref. [TT|). The 
AFM phases are the Neel phase (N) shown in Fig. [IJa) , 
the striped (S) phase discussed in our earlier paper [23], 
and the anti-Neel (aN) phase shown in Fig. QJc). The S, 
aN, and N states have, respectively, 1, 2, and all 3 NN 
spins to a given spin antiparallcl to it. Equivalently, if we 
consider the sites of the honeycomb lattice as comprising 
a set of parallel sawtooth (or zigzag) chains (in any one of 
the three equivalent directions), the S state comprises al- 
ternating FM chains, while the aN state comprises AFM 
chains in which NN spins on adjacent chains are paral- 
lel. Although there are infinite manifolds of non-coplanar 
states degenerate in energy with each of the S and aN 
states at T = 0, both thermal and quantum fluctuations 
[TT] select the collinear configurations. When J3 > there 
is a region in which the spiral state shown in Fig. [IJb) is 
the stable gs phase. It is characterized by a spiral angle 
defined so that as we move along the parallel sawtooth 
chains [drawn in the horizontal direction in Fig. QJb)] the 
spin angle increases by n + <p from one site to the next, 
and with NN spins on adjacent chains antiparallel. The 
classical gs energy is minimized for this spiral state when 
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We note that as <j> — > this spiral state becomes the 
collinear N state with energy per spin, 
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and there is a continuous phase transition between these 



two states on the boundary y 



r, for i < x < 



where y = J3/J1 and x = J2/ J\ ■ Similarly, as <fi —> tt the 
spiral state becomes the collinear S state, and there is a 
continuous phase transition between the two states on the 



boundary line y 



2 , for x > h . There is a first-order 



phase transition between the collinear N and S states along 
the boundary line x — |, for y > |. These three phases 
(N, S, and spiral) meet at the tricritical point (x, y) = 
(|, 5). As x — s- 00 (for fixed finite y), the spiral pitch angle 
4> — > |7r. In this limit the model becomes two HAFMs 
on weakly connected interpenetrating triangular lattices, 
with the usual classical ordering of NN spins oriented at 
an angle to each other on each sublattice. 
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The above three states are the only classical gs phases 
when y > 0. For y < the N state persists in a region 
bounded by the same boundary line as above, y = §#—3, 
for — 2 < x < g, on which it continuously meets a second 
spiral state, and by the boundary line y = — 1, for a; < 
— i, at which it undergoes a first-order transition to the 
FM state, which itself is the stable gs phase in the region 
x < — 5 and y < — 1. Another collinear AFM state, the 
aN state shown in Fig. [ljc), with energy per spin, 

= y(-Ji-2J 2 + 3J 3 ), (4) 

becomes the stable gs phase in the region x > ^, for y < 
\{x— [x 2 + 2(x~ I) 2 ] 1 / 2 }- On the boundary it undergoes a 
first-order transition to the spiral state shown in Fig. [TJb) . 
For i < x < i the spiral state shown in Fig. [IJb) meets 
another spiral gs phase on the boundary line y = 0, at 
which point there is a first-order transition. The pitch 
angle of this second spiral phase smoothly approaches the 
value zero along the above boundary with the N state, 
and the value n along a second boundary curve that joins 
the points (x,y) = (— |, — 1) and (5, 0), on which it meets 
the aN state. Both transitions are continuous ones. This 
second spiral meets the three collinear states N, aN, and 
FM at the tetracritical point (x,y) = (—5, — 1). 

Henceforth we restrict consideration to the J1-J2 model 
where J3 = (and J\ = 1). The classical (s — > 00) model 
thus has the N state as its gs for J 2 < gi whereas for 
J-2 > \ the gs comprises an infinite family of degenerate 
coplanar states with spiral order [including that shown in 
Fig. |Tfb)] , in which the spiral wave vector can point in 
any direction [TTJ[25j[26] . It is found [26] that, to leading 
order in 1/s, spin wave fluctuations lift this degeneracy 
by the well-known order-by-disorder mechanism, in favour 
of specific wave vectors. However, these spiral states for 
the J\-Ji model are expected to be very fragile against 
quantum fluctuations, and indeed to leading order in 1/s 
the spin- wave correction to the spiral order parameter has 
been shown to diverge as log TV, where N is the number of 
lattice sites [55], although it is still possible that higher- 
order terms in 1/s involving spin- wave interactions could 
stabilize the spiral order for large enough values of s. In 
view of the close proximity of the classical collinear AFM 
aN state for small values of the NNNN coupling J3 in the 
J1-J2-J3 model, it seems very likely that spiral order in 
the spin-i J1-J2 model might well be totally absent. 

To investigate this question further, and more generally 
to consider the entire T = phase diagram of the J\- 
J2 model, we utilize the coupled cluster method (CCM) 
[2"TH2"5] as in our work for the analogous J1-J2-J3 model 
with J 3 = J 2 [23] ■ When used at high orders in the sys- 
tematic approximation schemes developed for it, the CCM 
is an accurate approach to tackling a wide variety of quan- 
tum spin systems [5H71 I50H5B] . In particular, it can accu- 
rately locate the quantum critical points (QCPs) in such 
frustrated systems prni3"Tlf32,34,36]. as well as helping to 



determine the nature of the phases involved, including any 
quantum paramagnetic phases without magnetic order [7] . 

Since the CCM is a size-extensive method it provides 
results in the limit N — ¥ 00 from the outset. However, 
it requires us to input a model (or reference) state, with 
respect to which the quantum correlations may, in prin- 
ciple, be exactly included (and see, e.g., Refs. [51 l3"Tll3"8] 
and references cited therein). We use here the Neel (N), 
spiral, and anti-Neel (aN) states shown in Fig. [TJ as our 
CCM model states. The CCM then incorporates multi- 
spin correlations on top of the chosen gs model state $} 
for the correlation operators S and S that parametrize 
the exact gs ket and bra wave functions of the system 
in the respective exponentiated forms |\f') = e s |$) and 
(*| = ($|^e~ s , where (*|*) = 1. The Schrodinger kct 
and bra equations arc H\^) = E\^) and {^\H = £'(#| 
respectively. The correlation operators are defined as 
S = Y^i^a^i^t anc ^ S = 1 + "^Zi^qSiCJ respectively. 
The operators Cf = {Cj)^ and Cj are the creation and 
destruction operators respectively, where Cq = 1 and 
($|C/ = ; V7 ^ 0. The set {Cf = s+s+ • • • s+ }, where 
s^ = + is?, forms a complete set of multispin creation 
operators with respect to the model state |$) as a gener- 
alized vacuum. We then calculate the correlation coeffi- 
cients {<Sj,iSj} by minimizing the gs energy expectation 
value H = (^/\H\^f) with respect to each of them. This 
yields the coupled sets of equations (&\C]~e~ s He s \&) = 
and {<i>\S(e- s Hc s -E)Cf\<I>) = ; V/ ^ 0, which arc used 
to calculate the ket- and bra-state correlation coefficients 
within specific truncation schemes on the retained set {/} 
described below. It is necessary to use parallel computing 
routines for high-order computation [3TH55] . 

For the s = | case considered here we use the well-tested 
localized LSUBm truncation scheme which includes all 
multi-spin correlations in the CCM correlation operators 
over all regions on the lattice defined by m or fewer con- 
tiguous lattice sites. The numbers Nf of such fundamental 
configurations that are distinct under the symmetries of 
the lattice and the model state in various LSUBm approx- 
imations increases rapidly with the truncation index m. 
For example, the highest LSUBm level that we can reach, 
even with massive parallclization and the use of super- 
computing resources, is LSUB12, for which Nf = 293309 
for the aN state. The raw LSUBm data still need to be 
extrapolated to the exact m — > 00 limit. Although there 
are no exact extrapolation rules we have a great deal of 
experience in doing so. Thus, for the gs energy per spin, 
E/N, we use (see, e.g., Refs., [MM3I33I31I3E]) 

E(m)/N = ao + aim~ 2 + 02m -4 ; (5) 

while for the magnetic order parameter (sublatticc mag- 
netization), M, we use either the scheme 

M (m) = b + bxtrT 1 + b 2 m~ 2 , (6) 

for systems showing no or only slight frustration (see, e.g., 
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Fig. 2: Difference between the gs energies per spin (e = 
E/N) of the spiral and anti-Neel phases (Ae = e spiral - 
e aN ) versus Ji for the spin-i J1-J2 honeycomb model 
( Ji = 1) in LSUBm approximations with m = {6, 8, 10}. 



Fig. 3: CCM LSUBm results for the gs energy, E/N, of 
the Neel and anti-Neel phases of the spin-i J\-Ji hon- 
eycomb model (Ji = 1), with m = {6,8, 10, 12} and the 
extrapolated LSUBoo result using this data set. 



Refs. [3T)Il3T] ). or the scheme 

M(m) = c + cito~ 1/2 + b 2 m^ /2 , (7) 

for more strongly frustrated systems or ones showing a 
gs order-disorder transition (see, e.g., Refs. [6ll7l l34ll36] ). 
Since the hexagon is an important structural element of 
the honeycomb lattice we perform the extrapolations only 
for LSUBm data with m > 6. 

We show in Fig. [5] the difference in the gs energies of the 
CCM LSUBm results based on the spiral and aN model 
states. For the spiral state the pitch angle at a given 
LSUBm level is chosen to minimize the corresponding es- 
timate for the gs energy. Although the energy differences 
are small the results for all values of m, as well as the ex- 
trapolated results, show clearly that, as expected from our 
previous discussion, the spiral state that is the classical gs 
for all values n = J2/J1 > 0.5 gives way to the collinear 
aN state as the stable gs phase for the spin-i model out 
to much higher values of k. If a quantum phase transition 
between the spiral and aN states does exist, Fig. [5] shows 
that it can occur only at a value k > 1. Henceforth we 
concentrate on gs phases other than the spiral phase. 

In Fig. |3] we show results for the gs energy per spin, 
E/N, using the N and aN states as model states. We ob- 
serve that each of the energy curves based on a particular 
model state terminates at some critical value of k (that it- 
self depends on the LSUBm approximation used), beyond 
which no real CCM solution can be found. We note that in 
Fig. [U results are shown for each LSUB777, case down to val- 
ues of k at which real solutions based on the spiral model 
state cease to exist. In all cases the corresponding termi- 
nation point at a given LSUB7« level shown in Fig. [3] for 
aN model state is lower than that for the equivalent spiral 
model state case. Such terminations of the CCM solutions 
are well understood They are simply manifesta- 

tions of the quantum phase transitions in the real system, 
and may thus be used to estimate the positions of the 
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Fig. 4: CCM LSUBm results for the gs order parameter, 
M, of the Neel and anti-Neel phases of the spin-^ J\-Ji 
honeycomb model (Ji = 1), with m = {8, 10, 12}, and the 
extrapolated LSUBoo(l) and LSUBoo(2) results using this 
data set and Eqs. ^ and ([7]) respectively. 

corresponding QCPs [29], although we do not do so here 
since we have more accurate criteria available as discussed 
below. We note, however, that as is usually the case, the 
CCM LSUBm results for finite m values for both the N 
and aN phases shown in Fig. [3J extend beyond the corre- 
sponding LSUBoo transition points. For large values of m 
each LSUBm transition point is quite close to the actual 
QCP where that phase ends. For example, the LSUB12 
termination points shown in Fig. [3] are at nf ~ 0.23 for 
the N state and K aN 0.35 for the aN state. The CCM 
results show a clear intermediate regime in which neither 
of the quasiclassical AFM states (N or aN) is stable. 

We now discuss the magnetic order parameter, M , in or- 
der to investigate the stability of quasiclassical magnetic 
LRO. Our CCM results for M are shown in Fig.|U The ex- 
trapolated Neel order parameter goes to zero at a value k Ci 
that is very insensitive both to which extrapolation scheme 
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is used of Eqs. (JSJ) or (J7J, and to whether or not we include 
the LSUB6 data point in the extrapolations. Our best es- 
timate is K Cl « 0.207±0.003. This value can be considered 
as our first CCM estimate of the corresponding QCP of 
the model. It is in reasonable agreement with, but much 
more precise than, similar estimates for k Ci of about 0.17- 
0.22 from an exact diagonalization (ED) approach |16j . 
about 0.2 from an alternative ED approach [T3], about 0.2 
from a Schwinger boson approach jlOj . and about 0.13-0.17 
from a pseudo-fermion functional renormalization group 
approach [T5], but in substantial disagreement with a re- 
cent variational Monte Carlo (VMC) estimate of about 
0.08 [3D]. As the authors admit, the VMC study seems to 
substantially underestimate the QCP K Cl at which Neel 
order disappears. As expected, our own estimate shows 
that, as usual, quantum fluctuations preserve the collincar 
Nccl order to stronger frustrations than the corresponding 
classical transition to noncollincar spiral order at k c \ = g . 

By contrast with the situation at the lower QCP at k Ci , 
Fig. [4] shows that the corresponding QCP at k C2 at which 
the anti-Neel order vanishes is considerably more difficult 
to estimate from the extrapolated CCM LSUBm values, 
with estimates that range from 0.47 to 0.64. We find a 
much more accurate estimate for k C2 below. Neverthe- 
less it is clear already that a new quantum phase exists 
in the range k Ci < k < k C2 . It is also clear that, as 
suggested above, the two QCPs are very close to the cor- 
responding CCM termination points seen in Fig. [31 In our 
previous work on the spin-^ J1-J2-J3 model on the hon- 
eycomb lattice [53] with J\ = 1, we found strong evidence 
for a nonmagnetic plaquctte valence-bond crystal (PVBC) 
phase along the line J3 = J2 bewteen the two quasiclas- 
sical AFM phases, namely the Neel (N) and striped (S) 
phases. It seems likely that the corresponding phase in 
the present J3 = case, which now intervenes between the 
N and aN phases, might also be the same PVBC phase. 

In order to investigate the possibility of a PBVC phase 
we consider a generalized susceptibility \f that describes 
the response of the system to a "field" operator F (see, 
e.g., Ref. [7]). A field term F = 5 6 is added to the 
Hamiltonian (fTJ), where O is an operator which corre- 
sponds here to the possible PVBC order illustrated in 
Fig. [5] and which thus breaks the translational symmetry 
of H. The energy per site E(S)/N = e(5) is then calcu- 
lated in the CCM for the perturbed Hamiltonian H + F, 
using both the N and aN model states. The susceptibility 
is defined as xf = — (d 2 e(S))/(d5 2 )\ g=Q . Clearly, the gs 
phase becomes unstable against the perturbation F when 
Xp 1 becomes zero. As in Ref. [53] we use the extrapolation 
scheme X^ 1 (m) = do + d\mT 2 + d2m~ 4 . 

Our CCM results for Xp 1 are shown in Fig.[S] The num- 
ber of LSUB12 fundamental configurations for the pla- 
quctte susceptibility for the aN state is Nf = 877315. 
The extrapolated inverse susceptibility vanishes on the 
Neel side at K « 0.24 ± 0.01 and on the anti-Neel side 
at k ~ 0.385 ±0.010. The shape of the CCM curves where 
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Fig. 5: (Color online) Left: CCM LSUBjti results for the 
inverse plaquette susceptibility, 1/x, of the Neel and anti- 
Neel phases of the spin-i J\-Ji honeycomb model (Ji = 
1) with m — {6,8,10,12}, and the extrapolated results 
LSUBco(l) using m = {6, 8, 10, 12} and LSUBoo(2) using 
m = {8,10,12} (see text). Right: The fields F = 5 6 
for the plaquette susceptibility x- Thick (red) and thin 
(black) lines correspond respectively to strengthened and 
weakened NN exchange couplings, where O = Y^u j\ • 
Sj, and the sum runs over all NN bonds, with ay = +1 
and —1 for thick (red) and thin (black) lines respectively. 
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on the Neel side is strongly suggestive of a contin- 



uous transition there, just as we found for the correspond- 
ing J1-J2-J3 model [53]. The shallow slope of the Xi? 1 
curves there makes it difficult to estimate accurately the 
point where it vanishes. Nevertheless it is certainly con- 
sistent with the much more accurate value we obtained for 
k Ci above from the point where the Neel order parameter 
M vanishes. On the other hand we cannot exclude the 
possibility of the transition between the Neel and PVBC 
states occurring via an intermediate phase that exists in 
the region K C1 < k < 0.24. Just such an intermediate (res- 
onating valence bond spin-liquid) state has been discussed 
in Ref. [11]. By contrast, the shape of the CCM curves 
for Xp 1 on the anti-Neel side are much more indicative of 
a first-order transition, and the point where Xf 1 ~~ ^ ^ on 
that side gives us our best estimate for k C2 w 0.385±0.010. 
Again, this value is in good agreement with, but much 
more accurate than, estimates of about 0.35-0.4 from two 
different ED studies [T31[To1 . On both the N and aN sides it 
seems very likely that the PVBC phase occurs at, or very 
close to, the QCPs where the quasiclassical magnetic LRO 
in the N and aN phases vanishes. Since the N and PVBC 
phases break different symmetries, and our CCM results 
show that they appear to meet at n Cl w 0.21 at a contin- 
uous transition, they support the deconfinement scenario 
there. The possibility of deconfined quantum criticality 
for the frustrated honeycomb HAFM was pointed out in 
Ref. [T], and supporting evidence for the J1-J2-J3 model 
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has also been found both by us [53] and by others [15] . 

We have studied the influence of quantum spin fluctu- 
ations on the gs properties of the spin-i J\~Ji HAFM 
(Ji > Q,K = J2/J1) on the honeycomb lattice. We find a 
paramagnetic PVBC phase in the regime k Ci < k < k c . 2 , 



where n r 



0.207 ± 0.003 and k c 



0.385 ±0.010. The 



transition at k Ci to the Neel phase appears to be a con- 
tinuous transition of the dcconfincment variety, while that 
at k C2 is of first-order type to an anti-Neel-ordered AFM 
phase that does not occur in the classical version of the 
model except at the critical point k = ^. Our results in- 
dicate that the spiral phases that exist classically for all 
values k > i are absent for all values k < 1 , but may exist 
for larger values. To investigate this and other aspects of 
the model further we shall present results in a future pa- 
per on the phase diagram of the extended J1-J2-J3 model, 
using the same CCM techniques. 

We thank the University of Minnesota Supercomputing 
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